Quantum many-body simulations using Gaussian phase-space representations 
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Phase-space representations are of increasing importance as a viable and successful means to study exponen- 
tially complex quantum many-body systems from first principles. This paper traces the background of these 
methods, starting from the early work of Wigner, Glauber and Sudarshan. We focus on modern phase-space 
approaches using non-classical phase-space representations. These lead to the Gaussian representation, which 
unifies bosonic and fermionic phase-space. Examples treated include quantum solitons in optical fibers, collid- 
ing Bose-Einstein condensates, and strongly correlated fermions on lattices. 
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I. INTRODUCTION 

In this paper, we will trace how the concept of coherence 
and coherent states has led an important advance: the quan- 
tum phase-space representation. Through the development of 
phase-space representations, the idea of coherence can help 
the understanding and simulation of the physics of many-body 
systems, both in thermal equilibrium, and in time-dependent, 
quantum dynamical calculations. This is of increasing impor- 
tance beyond quantum optics, as new experiments explore the 
quantum correlations and dynamics of interacting particles. 

We show that a more general approach to coherence leads 
to the Gaussian phase-space method, which unifies the repre- 
sentation of both bosonic and fermionic many-body systems. 
This powerful idea has many ramifications. It encompasses 
all the known bosonic representations in a simple, clear for- 
malism, and extends these ideas to fermions as well. It is also 
extremely useful in applications, as we will show using both 
equilibrium and non-equilibrium examples. 

A particular quantum state that illustrates this is the coher- 
ent state. It was introduced originally by Schroedingerfjl for 
the harmonic oscillator, and later applied to the radiation field 
through the seminal work of Sudarshan|2] and Glauber|3]. 
These states are fully coherent in the sense that normally or- 
dered operator moments factorize to all orders. 

The definition of a coherent state is extremely simple. If a 
is a field-mode annihilation operator, then the coherent state 
is defined as a normalized eigenstate of a, 



Husimi|5] Q-function. The closely related operator associ- 
ations of Lax 1 6], Agarwal|0 @] and co-workers were used 
to develop a quantum theory of the laser. While useful for 
the laser, these all lack essential ingredients that would al- 
low them to be useful as a probability distributions in first- 
principles many-body dynamical simulations. Most are sim- 
ply non-positive, as in the case of the P-function and Wigner 
function. Any representation that uses a classical-like phase 
space has no corresponding exact stochastic equation when 
there are inter-particle interactions. 

We will explain how this problem is solved by ex- 
tending the phase-space dimension, giving rise to the 
positive P-representationOJ [Toll . A unifying principle 
is the use of non-orthogonal basis sets, which leads to 
the idea of a stochastic gauge svmme try |l l |. and more 
general Gaussian phase-space methods |12J uM- These 
have many applications to interacting Bose and Fermi 
systems. Both thermal equilibrium and first-principles 
quantum dynamical time-evolution (either unitary or dis- 
sipative) can be treated. Recent bosonic examples include 
quantitatively tested predictions on quantum soliton time- 
evolution lll4ll . as well as novel predictions for topical exper- 
iments including: colliding Bose-Einstein condensates| 15], 
tunnel-coupled condensates||Tq|. sup erchemistry fl7ll . 
molecular dissociation! 1 4 Il9l l20t l2lll . micro-mechanical 
resonators 1 23 , triple EPR correlations] 23, and non- 
equilibrium criticality in parametric downconversion l24i . We 
also give results for phase-space simulations of the fermionic 
Hubbard model in thermal equilibrium! 25]. 
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These states form a complete mathematical basis, providing 
examples of quantum states which are perfectly coherent to 
all orders. The idea can be extended to other algebras, for 
example the SU(N) coherent states, and were used to construct 
the P-representation - a representation of the radiation field in 
terms of diagonal coherent state projection operators. This 
quantum operator representation has the form (for a single 
mode) of: 



P(a) \a) (a 



d 2 a. 
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This representation maps a quantum state into a distribu- 
tion on a classical phase space. Other representations like 
this exist, including the Wignerfl representation and the 



II. QUANTUM MANY-BODY SYSTEMS 

Quantum many-body theory is the generic theory we cur- 
rently use for describing all non-astronomical physical sys- 
tems from a microscopic point of view. It is applicable to a 
wide range of problems. 



A. Ultra-cold atomic Gases 

As simple examples of interacting quantum systems, con- 
sider the ultra-cold atomic Bose-Einstein condensates and de- 
generate atomic Fermi gases. Ultra-cold atoms are an ideal 
quantum many-body system. In these experiments, the in- 
teracting atoms are isolated from other matter, by virtue of 



2 



being optically or magnetically trapped in a high-vacuum en- 
vironment at low temperatures. Important advances in the 
last decade include: Bose-Einstein condensation (BEC), atom 
lasers, superfluid Fermi atoms, superchemistry (stimulated 
molecule formation), atomic diffraction, interferometers, and 
temperatures below InK. 

Such well-controlled and simple physical systems present 
an opportunity to quantitatively test quantum mechanics in 
new regimes, where macroscopic and many-body effects play 
a dominant role. 



B. Many-body quantum dynamics 



Exact solutions - even if all the energy eigenstates are known 
(which is unusual) evaluating the initial expansion co- 
efficients for quantum dynamics remains exponentially 
difficult, and therefore impractical 

New hardware - Feynman proposed quantum computers to 
solve many-body problems - currently, these do not ex- 
ist beyond 2 — 4 qubit capacity 

New software - Gaussian quantum phase-space simulation 
methods can give practical techniques using existing 
computers, simulating quantum systems equivalent to 
nearly a million qubits. 



Before one can make quantitative predictions, there is a sig- 
nificant problem to overcome: quantum many-body problems 
are exponentially complex. 

To illustrate this, consider a Bose gas with N atoms dis- 
tributed among M modes. Each mode can have one or all 
atoms. The number N s of quantum states available is: 



N» = 



(n + m -iy. 

N\ (M- 1)! 



(3) 



A typical BEC may have N 
astronomical number of: 



M 



500, 000, giving the 



N, = 2 2N = 10 



-,300,000 



(4) 



Hilbert space dimension can also be classified by the num- 
ber of equivalent quantum bits (qubits), which is log 2 N s — 
2N = 1, 000, 000, in this example. 

There are a number of possible solutions to dynamical prob- 
lems. Here we focus on methods which are exact, in the sense 
that errors can be estimated and reduced where necessary. 
As an example, while Density Matrix Renormalisation Group 
(DMRG) methods [26] can be useful for one-dimensional cal- 
culations, including dynamics, the Hilbert-space truncation 
is not always a well-controlled approximation. Similar dif- 
ficulties occur in the density functional approach ll27ll . Uncon- 
trolled approximations cannot be used as a basis for testing 
quantum mechanics. Any discrepancies observed may sim- 
ply be caused by calculational errors, rather than fundamental 
issues. 

Candidates for exact solutions are as follows: 

Path integrals and Monte-Carlo - these are useful for 
bosons at thermal equilibrium. For quantum dynam- 
ics and for fermions, there are phase and sign problems, 
making these methods often impractical. 

Perturbation theory - while applicable for certain problems, 
this method generally doesn't converge in quantum field 
theory 

Numerical diagonalization - the problem of an exponen- 
tially large matrix size rules out such brute force meth- 
ods, except for very small particle numbers 



III. QUANTUM PHASE-SPACE METHODS 

The great power of phase-space methods is their ability to 
accurately compute the quantum dynamics of fully macro- 
scopic systems directly from the Hamiltonian, without resort- 
ing to overarching approximations. This confers several ad- 
vantages over previous methods, despite the introduction of 
randomness that limits precision: 

Firstly, all uncertainty in the results is confined to random 
statistical fluctuations, with no systematic bias. Impor- 
tantly, the magnitude of this uncertainty can be reliably 
estimated from the distribution of sub-ensemble means 
by using the central limit theorem 

Secondly, these methods lead to relatively simple equations 
that can be easily adapted to trap potentials and local 
losses, whose magnitude and shape can be chosen arbi- 
trarily. This is in stark contrast to approximate methods, 
which can become much more complicated or even in- 
applicable under such conditions. 

The Gaussian quantum phase-space representation described 
here encompass all the earlier known phase-space methods. 
Therefore, we start by reviewing these earlier approaches. 



A. Classical and Quantum phase space 

Wigner [4] originated the idea of a classical-like phase- 
space or quasi-probability description. For M modes, these 
methods scale linearly with mode number, having just M 
complex dimensions. Variations on this theme include the 
Husimi Q-function|5], and the Glauber-Sudarshan|2j |3| P- 
representation. For many quantum states, they result in a 
positive-valued distribution. For the Q-function, this is always 
true. Despite this, one finds that there is no corresponding 
stochastic equation for cubic or quartic Hamiltonians. Thus, 
there is no method for efficiently time-evolving a sampled dis- 
tribution of an interacting system, except through an approxi- 
mate truncation of the equations of motion. 

The solution to this problem is to use an enlarged phase 
space, which includes off-diagonal terms in a coherent-state 
expansion. Intuitively, this allows for quantum superpositions 
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Figure 1 : The full variance a p is composed of a distribution variance 
dp, and a basis variance <ta- 

between more than one classical configuration. The simplest 
possibility is the positive-P (+P) distribution! 10], which has 
2M coordinates. It results in a distribution function which is 
always positive, and given certain conditions, obeys a stochas- 
tic equation. It has the definition that: 

p= f P( ai (3)\^\d 2 ad 2 f3. (5) 
J (/311a) 

B. Quantum phase-space representations 



and Glauber showed, can be used to defined phase-space rep- 
resentations for fermions|31]. Unlike their bosonic coun- 
terparts, the fermionic coherent states have no direct phys- 
ical meaning. Moreover, while they are useful for formal 
calculations, they have limited applicability as a basis for 
practical, numerical calculations, because the of the com- 
plexity that arises from the anticommuting properties of the 
algebra[32|. The coherent-state P-representation is then a 
function of Grassmann numbers, not a probability. 

But this anticommuting complexity is related to the unphys- 
ical states contained in the coherent basis. Fermionic coherent 
states require Grassmann numbers because of the way they 
include coherences between states with an odd number dif- 
ference. Consider a coherent superposition of zero and one- 
particle states:]^) = a |0) + /3|1), which gives a nonzero 
value for the coherence (a) = a* j3. Because the |1) state in- 
volves an anticommuting operator, one of the amplitudes must 
also be anticommuting, for consistency. This also means that 
the coherent amplitude is anticommuting. 

However, from superselection rules, we know that fermions 
can only be created in pairs, and thus such superpositions are 
excluded. Thus one can avoid this anticommuting problem by 
considering an operator basis which only includes coherences 
that are allowed by the superselection rule. 



Guided by the formalism of Equations (|2j and 0, one can 
define a general quantum phase-space representation by ex- 
panding the density matrix 2 using a complete basis of opera- 
tors A( A ): 



2 = / P( A )A( X)d\ 



(6) 



Provided P{ A ) remains positive and sufficiently bounded, 

quantum dynamics can be transformed into trajectories in A . 

Different basis choices for A( A ) then result in different rep- 
resentations. For example, the P-representation has a single 
complex dimension (for M = 1), so Ai = a, and: 



A(a) = \a) (a\ . 



(7) 



As shown in Figure Q, there are trade-offs in the choice 
of basis, since the quantum variance is partly due to the dis- 
tribution, and partly due to the basis. By minimizing the the 
distribution variance, one can reduce the sampling error of the 
representation. This typically involves an over-complete, non- 
orthogonal basis in which each member of the basis is closely 
matched to a physical state that occurs in the simulation. 



C. Fermionic phase space 



D. General M-mode Gaussian operator 

The most general phase-space representation, for both 
fermions and bosons, is obtained with Gaussian operators. 
These provide an (over)complete basis for fermions even 
when the coherences, and thus Grassmann components are 
excluded| 33]. These also generalize the concept of coherence: 
physical states with Gaussian density operators have operator 
products that factorize in a similar, but more general way than 
coherent states. 

To define these, we introduce a as a column vector of M 
bosonic/fermionic annihilation operators (indicated as the up- 
per or lower sign respectively), and 2^ the corresponding row 
vector of creation operators,. Their commutation relations are: 



a-k, a 



(8) 



A Gaussian operator is defined as a normally ordered ex- 
ponential of a quadratic form in annihilation and creation 
operators. Introducing extended 2 M- vectors of operators: 
2 = (2, (2t) T ), with adjoint defined as a} = (2t,2 T ), the 
operator fluctuation is then: Sa = 2 — a , where a = (a, (3*) 
is a 2M-vector c-number. A Gaussian operator can therefore 
be written as: 



Coherent states for fermions|28, 29] can be defined by 
means of anti-commuting Grassmann numbers, and have 
been used, for example, in path-integral formulations for 
fermions l30ll . Like their bosonic counterparts, fermionic co- 
herent states provide an overcomplete basis set, and as Cahill 



A, 



fi |cx| =Fl ^ 2 : exp 



<52 f | / 



=/ 2= 



(Si I 



■ (9) 



In the fermionic case the square root of the determinant (for 
normalization purposes) is to be interpreted as the Pfaffian of 
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the matrix, in an explicitly antisymmetric form. The addi- 
tional factor J in the exponent only appears in the fermionic 
case: 



±1 
I 



E. Operator mappings 



(10) 



The covariance a is best thought of as a kind of dynamical 
Green's function. It can be expanded as: 



±n T 
m+ 



m 
n 



(11) 



Here n is a complex matrix whose average is the normal 
Green's function for particles, while n = l±n. In many-body 

terminology, m and m + correspond to anomalous Green's 

— > 

functions. The representation phase space is therefore A = 
(f2, ct, (3, n, m, m + ) for bosons; in the case of fermions, one 
must set a = (3 = 0. 

The significance of the definition of n and m is that it leads 
to useful bosonic and fermionic operator identities. For exam- 
ple, one finds that: 



= + n ij) P , 

where the weighted average is defined as: 



(12) 



O 



0(X 



o(\)np(x,T)d\ 



(13) 



For representations with fixed riy, one thus obtains a gener- 
alized operator-ordering. Classical phase-space distributions 
are recovered on setting on — (3* , and riy = cSij . For exam- 
ple, the Glauber-Sudarshan P-representation has c = 0, while 
the Wigner distribution has c = —1/2. More generally, this 
type of phase space allows for a stochastic covariance, which 
can dynamically change in time and space to suit the physical 
system. 

Other useful identities involve the relationship between the 
action of operators on the kernel, and the corresponding dif- 
ferential operators acting on the distribution itself. For sim- 
plicity, these are given in the number-conserving case (a = 



(3 = 0, m = 0, ra H 

nA 
An 



where 



(?/fln) 



0): 



nP-(I±n)— nP (14) 
an 

nP - n— (I ± n) P . 
on 



d /driji is a differential operator that 



F. Evolution equations 

There are three main types of problems studied with this 
approach, which provides a unified method for interacting 
fermions and bosons: 

• Canonical ensembles - thermal initial conditions 

• Quantum dynamics - unitary nonlinear time-evolution 

• Master equations - open system time-evolution to a 
steady-state. 

The purpose of the phase-space representation is to trans- 
form exponentially complex operator equations into tractable 
phase-space equations, which can then be effectively sampled 
via probabilistic means. For example, suppose that we wish 
to calculate a thermal ensemble. The grand-canonical density 
operator can be written as an operator differential equation, 



dp 



H - p,N , p 



L[P] 



(15) 



Similarly, one can also treat unitary evolution or evolution 
under a master equation as a generalized Liouville operator. 
By making use of the operator identities above, and provided 
conditions of compactness that allow partial integration are 
satisfied, one can transform the exponentially large operator 
equation into a stochastic equations that can be treated either 
numerically or, in some cases, even analytically. The generic 
form that results, in the Ito calculus, is: 



dn/dt 

dX/dt 



Q[U+g 
A + B(C- 



g), 



(16) 



acts both to the left and the right. 



where £ is a vector of Gaussian white noises. The function g 
is a 'stochastic gauge' function, that can be adjusted to guar- 
antee the stability of the resulting drift equations. 

In summary, this method greatly extends the approaches of 
Glauber, Sudarshan, Husimi and Wigner. No approximations 
are needed, apart from the sampling error, which can be es- 
timated and reduced by using more samples. The represen- 
tations use positive, nonsingular distributions on a relatively 
small (non-exponential) phase space. This reduces the overall 
complexity enormously. The price that is paid is that many 
trajectories can be needed to control sampling error, which 
typically grows with time. One must also design an appro- 
priate stabilizing gauge g, as stable trajectories are essential 
to remove boundary terms. The overall procedure is outlined 
schematically in Fig|2] 



IV. BOSONS 

The simplest general model of an interacting Bose gas is 
the Bose-Hubbard model, which includes nonlinear interac- 
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Figure 2: Strategies that need to be considered and optimized in 
quantum simulations. 



tions at each site, together with linear interactions coupling 
different sites: 



H{a, a+) = h Wi 3 -oJ°j + 



(17) 



where the frequency termwy is a nonlocal coupling, which 
includes chemical potential. The boson number operator is 
rii = o|tij . The most commonly used technique here is 
the positive-P representation, although more general Gaussian 
methods are also possible. 



A. Single-mode phase-diffusion 

As an example, consider the case of a single potential well 
containing a BEC in an initial coherent state. After apply- 
ing the relevant operator mappings, one obtains the following 
time-evolution equations: 



. da 

.dp 
~dr 

dn 



Re [(3a}+oj + v^~i( 2 (T) 

tt[<?lCl(T) + <?2C2(V)] ■ 



a 







(18) 



Here, unitary evolution leads to nonlinear phase-diffusion, 
as has been experimentally observed|34|. The stochastic tech- 
nique can be utilized to carry out a simulation of quantum 
evolution of an initial coherent state of up to 10 23 bosons! 

This is shown in Fig. [3] where first 100 atoms, and then 10 23 
atoms were simulated after appropriate choices of gauges gi% 
and noises d.2 1 35 ] . A time-reversal test of unitary evolution 
was carried out by reversing the sign of the Hamiltonian, in 
order to observe a recurrence to the initial physical state. This 
is even possible experimentally, using Feshbach resonances to 
control the interaction. 

The distribution graphs in Fig.|4]demonstrate that the mech- 
anism for the recurrence is not through a recurrence of the en- 
tire distribution, as only the physically observable moments 
have to show recurrence. The non-uniqueness of the basis 




x 10 




Figure 3: Simulation of (a) 100 and (b) 10 23 atoms in a single-mode 
trap, showing phase-decay together with a recurrence due to time- 
reversal. 



means that the final distribution is actually different to the ini- 
tial one; the effect of time-reversal is to change the detailed 
structure of the diffusive broadening, so that the final and ini- 
tial distributions have an equivalent physical density matrix. 



B. Optical fibre squeezing experiment 

To a very good approximation, photons in an optical fibre, 
with the Kerr nonlinearity present, are an experimental im- 
plementation of the famous one-dimensional Bose gas model 
in quantum field theory 113a. 13711 . Phase-space methods were 
used to make first-principles, testable predictions of quantum 
squeezing in this environment. We will show that these results 
are in excellent quantitative agreement with experiment, even 
including dissipation. 

We focus on recent polarisation squeezing experiments! 38], 
which are an efficient and flexible method for generating 
quantum states in the fibre| 14]. The experimental set-up is 
illustrated in Figure[5] Pulses are generated in pairs and prop- 
agate down orthogonal polarisation modes of an optical fibre. 
They are then combined in a Stokes measurement of polari- 
sation squeezing by means of a polarisation rotator, a beam 
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Figure 4: Phase-space distribution in the single-mode trap simula- 
tions with TV = 100 showing (a) time-reversal, (b) no time-reversal. 
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Figure 5: Polarisation-squeezing experiment! 36 



splitter and two detectors. 

Because the experiment involves ultrashort pulses, the 
quantum description must use photon-density operators 
V? x (t, z) and fy y (t, z) that include a range of spectral com- 
ponents: 



1 



dka a (t,k)e lik - ko)z+luJot , (19) 



where a — x,y. The commutation relations of these operators 
are $ a (t,z),$l,(t,z') = 6(z - z')5 aa , . 

For convenience, we use scaled variables in propagative 
frame: r = (t — z/v), ( = z/zo and (j> a = fy a y/vto/ri, 
where to is the pulse duration, zo = t§/|fc"| is the dispersion 
length and 2n = 2\k"\Ac/ {nzfkJ^to) is the photon number in 
a soliton pulse. 



To describe the evolution of the photon flux <j) a {T, £), we 
employ a quantum model of a radiation field propagating 
along a silica fibre, including x' 3 ^ nonlinear responses of the 
material and non-resonant coupling to phonons|3jl l4(ill .The 
phonons provide a non-Markovian reservoir that generates 
additional, delayed nonlinearity, as well as spontaneous and 
thermal noise. Because of fibre birefringence, the two polari- 
sation components do not temporally overlap for most of the 
fibre length, and so the cross-polarisation component of the 
Raman gain is neglected. The result, after discretization, is a 
Hubbard model like Ea.( ll7> . except with additional coupling 
to phonon reservoirs. 

The quantum operator equations are obtained by integration 
of the Heisenberg equations for the phonon operators to derive 
quantum Langevin equations for the photon-flux field: 

q^MtX) = |^2^(r,C)+if CT (r,C)^(T,C) (20) 

/oo 
dr'h(T - r')^(r', QMt', C)Mt, 0- 
-oo 

where the nonlinear response function h(r) includes con- 
tributions from both the instantaneous electronic response 
and the Raman response determined by the gain function 
The correlations of the reservoir fields 



a R {uj) 
are: 



(ft^cOf^c)) = ^^i[ nth (| w |) + e(-c)] 

X*(C-0*(w-W%a' ,(2D 

where is the temperature-dependent Bose distribution of 
phonon occupations. The Stokes (uj < 0) and anti-Stokes 
(u> > 0) contributions to the Raman noise are included by 
means of the Heaviside step function O. 

In all, we have over 10 8 photons in more than 10 2 modes, 
corresponding to an enormously large Hilbert space. Quan- 
tum dynamical simulations of such systems have been per- 
formed exactly using the +P representation 13(1 l43ll . How- 
ever, for large photon number n and short propagation dis- 
tance L, these exact squeezing predictions agree with a trun- 
cated Wigner phase-space method|44], which allows faster 
calculations. In effect, the Wigner representation maps a field 

operator to a stochastic field: f/) CT (C, T ) ~^ ^(C, r )- Stochas- 
tic averages involving this field then correspond to symmet- 
rically ordered correlations of the quantum system. Because 
of the symmetric-ordering correspondence, quantum effects 
enter via vacuum noise. 

After the mapping, we obtain a Raman-modified stochastic 
nonlinear Schroedinger equation for the photon flux that is of 
exactly the same form as Eq. ( 1201 139l 14511 . The correlations 
of the Raman noise fields T a and the initial vacuum noise are, 
respectively, 



(r CT ( w ,c)rv(u/,c')) = 



x<y(c-cO*(«-w% 



(Ai(r,0)Ae(r',0)) = — <5(r-r')<W. (22) 
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Figure 6: Antisqueezing and squeezing for L\ = 13.4m (squares) 
and Z/2 = 30m (diamonds) fibres. Solid and dashed lines show the 
simulation results for L\ = 13.4m and L2 = 30m, respectively. 
Dotted lines indicate sampling error in simulation results. Simula- 
tions are adjusted for linear loss of 24% and low-frequency GAWBS 
noise, which mainly affects the squeezing only at low power. Pa- 
rameters are parameters: to = 74fs, 20 = 0.52m, n = 2 x 10 8 , 
E„ = 54pJ and Ao = 1.51/Mn. 



Because of the symmetrically ordered mapping, the Stokes 
and anti-Stokes contributions to the Wigner Raman noise are 
identical. 

Antisqueezing and squeezing results are shown figure[6] for 
13.4m and 30m of fibre, with and without the excess phase 
noise included. The theoretical results for both squeezing and 
antisqueezing closely match the experimental data. The re- 
sults also show a deterioration of squeezing at higher intensity 
due to Raman effects, especially for longer fibre lengths. 



C. BEC collision with 150,000 atoms from first principles 

The collision of pure 23 Na BECs, as in a recent experiment 
at MIT|46], represents another opportunity for observational 
tests of first-principles quantum dynamical simulations lll5l 
HTn . In the simulations, a 1.5 x 10 6 atom condensate is pre- 
pared in a cigar-shaped magnetic trap with frequencies 20 Hz 
axially in the "X" direction, and 80 Hz radially ("Y" and "Z"). 
A brief Bragg laser pulse coherently imparts an X velocity of 
2vq w 20 mm/s to half of the atoms, which is much greater 



than the sound velocity of 3.1 mm/s. Another much weaker 
pulse generates a small 2% "seed" wavepacket at a Y velocity 
of v s = 9.37 mm/s relative to the center of mass. 

At this point the trap is turned off so that the wavepack- 
ets collide freely. In a center-of-mass frame, atoms are scat- 
tered preferentially into a spherical shell in momentum space 
with mean velocities v s s vq. Simultaneously, a four-wave 
mixing process generates a new coherent wavepacket at Y 
velocity -v s , as well as growing the strength of both of the 
wavepackets at ±v s by Bose enhanced scattering. 

The dynamics of the distribution of atom velocities and cor- 
relations between the scattered atoms have been calculated 
and are shown in Figures [7\ and [3] Such correlations have 
recently become experimentally measurable EH 14^. l50ll . and 
correlation behaviour qualitatively similar to that predicted by 
this model have been seen|51]. The system is described by: 



H 



2m 



V* f V$ 



^$t2$2 

2 



(23) 



The simulation is carried out using the positive-P 
representation in the center-of-mass frame from the 
moment the lasers and trap are turned off (t = 0). 
The initial wavefunction is modeled as the coherent- 
state mean-field Gross-Pitaevskii (GP) solution of the 
trapped t < condensate, but modulated with a factor 
[V0A9e tmv ^ h + Vfl49e- imv « x /' i + VoMe- imv <> y / H ] 
which imparts the initial velocities. The field Hamiltonian 
is discretized with a lattice size of 432 x 105 x 50, again 
generating a Hubbard-type Hamiltonian like Eq ( fTTt . 

As one might expect of a method that attacks such an expo- 
nentially complex problem, there are limitations. Most sig- 
nificantly, the size of the sampling uncertainty grows with 
time, and eventually reaches a size where it is no longer prac- 
tical to produce enough trajectories to get useful precision. In 
the above case a useful observable-to-noise ratio lasted un- 
til / < 410 [is. In general the simulation time possible de- 
pends on several factors: coarser lattices, weaker interactions, 
or smaller density all extend it. This time can be estimated 
using the formulae found in 15211 . Comparisons were made 
with a previous approximate simulation|53], using a truncated 
Wigner method 14 14411 . The approximate method was less ac- 
curate at large momentum cutoff, due to a diverging truncation 
error. 

The model treats M — 2.268 x 10 6 interacting momentum 

modes. Since each of the N = 1.5 x 10 5 atoms can be in 

any one of the modes, the Hilbert space contains about N s s» 

M N 10 1 ' 000 ' 000 orthogonal quantum states. In terms of 

accessible states at fixed number, there are w (M/N) s» 
2 60o,ooo stateS) or 60Qj qqq qubits 

This is the largest Hilbert space ever treated in a first- 
principles quantum dynamical simulation. 
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Figure 7: Momentum space snap-shots in the center-of mass frame, 
left: Velocity distributions in the axial (x) and radial (y) direc- 
tions, right: Radial distribution at x=0. The formation of the 
fourth coherent wavepacket at v y « — v s and the scattered shell at 
\v\ ~ v s = 9.37mm/s are seen. Logarithmic color scale. Average 
of 1492 trajectories. 



V. FERMIONS 

A. Hubbard model 

To demonstrate the utility of this fermionic representation, 
we next consider the fermionic Hubbard model 1 54]. This 
is well-known in condensed matter physics as the simplest 
model of interacting fermions on a lattice: 

H(n u n^) Ujn ijtlT + UJ2'. n 3ijA n, hj . x : . (24) 



Here <, = a\ a aj^, for lattice index j and spin index a = 
(|, J.) = (— 1, 1), while tij is the inter-site coupling and U is 
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Figure 8: Correlations between scattered atoms: time evolu- 
tion. Plate a shows the extremely strong number correlations 
g (vo, —vq) between atoms with opposite velocity (solid line) in 
the scattered shell at \vq\ = v s (away from the coherent wavepack- 
ets), and thermal correlations g^(vo, Do) = 2 between scattered 
atoms at the same velocity (dashed). Triple lines indicate uncertainty. 
Plate b shows the coherence width in velocity space for scattered 
atoms at similar velocities centered around vq. Plotted is the Full- 
width at half-maximum (FWHM) of (vo, v a + v)\. 
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Figure 9: ID Hubbard model: Correlation function (n^nf) versus 
temperature. The 100-site numerical solution is compared with the 
zero-temperature exact solution of an infinite lattice|60]: t — 1 and 
(7 = 2. 



the strength of on-site interaction between particles. 

Thought to be relevant to high-T c superconductors! 55j, the 
Hubbard model has had renewed interest because it describes 
an ultra-cold gas in an optical lattice|56], as has been experi- 
mentally realised by Kohl et al|57]. Within this simple model 
is a great complexity, that leads to sampling error problems 
for quantum Monte Carlo (QMC) methods because of nega- 
tive weights|58|. Such sign problems occur for repulsive in- 
teractions away from half-filling in two or more dimensions, 
and increase with lattice size and interaction strength|59]. 

Results of recent phase-space numerical simulations |25J in 
one and two dimensions are shown in Figures [9] and ^| The 
sampling error remains well-controlled at low temperatures, 
even for filling factors in 2D for which other QMC methods 
suffer sign problems. 

To explain the method used, we first note that the Hubbard 
Hamiltonian conserves number, so the number-conserving 
subset of Gaussian operators provides a complete basis, i.e. 
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Figure 10: 2D Hubbard model on a 16x 16 lattice: Energy as a func- 
tion of of temperature for various chemical potentials, t = 1 and 

[7 = 4. 



Associated with each stochastic path is a weight, governed by 
dVl/dr = —VlH (ni, n_i). Importantly, because the choice of 
mapping, the phase-space equations are real and the weights 
thus remain positive, avoiding the usual manifestation of the 
sign problem. 

More precise numerical simulations by Assaad et al|61] 
have revealed that there is difficulty in sampling ground state 
properties with these phase-space equations. However, they 
also show that the correct ground-state results can be obtained 
by a supplementing the phase-space simulations with a sym- 
metry projection procedure. 

Finally, we remark that the mapping from the Hubbard 
model to phase-space equation is far from unique. Thus these 
phase-space simulations of the Hubbard model may well be 
improved by appropriate choice of basis subset and stochastic 
gauge, as for bosonic simulations. 



the anomalous variables remain zero. The mappings given 
above can be applied to the grand canonical equilibrium equa- 
tion to give the Ito phase-space equations fPjll : 

^ = i{(I-n CT )Ti 1 W + n CT T( 2 )(I-„ (T )}.(25) 
The propagation matrices are defined for U > 0, as 

T&l = Ui ~ ki { Unjj,^ — fi + atV } , (26) 

where the stochastic terms are Gaussian white noises with the 
correlations 

(fj p) (r)#V)) = 2U5{T-r')8 j , j/ 6 r , r , . (27) 



VI. SUMMARY 

In summary, coherence theory and coherent-state methods 
leads to a unified phase-space representation for bosonic and 
fermionic quantum many-body systems, which are useful in 
simulations both in real time and in inverse temperature. Cal- 
culations have been carried out in one, two and three dimen- 
sions, with up to 10 23 particles and 10 6 modes. This is equiva- 
lent to a Hilbert space of nearly a million qubits. Phase-space 
ideas are also applicable to other complex sy stems 1 64 163[ - 
ranging from genetics, astrophysics, and biochemistry, to con- 
densed matter, particle physics and possibly even molecular 
physics. 
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